library(tidyverse)
library(metaumbrella)
collapsunique <- function(x) paste(unique(sort(x)), collapse = ", ")

source("1b - cct_level_cam_scoring_140424.R")
source("2 - cct_level_homogeneization_140424.R")
source("3 - ma_level_scoring140424.R")

chemin = "C:/Users/Corentin Gosling/drive_gmail/Recherche/Article 2 - Base de Donnees/7 - Data analysis/data/"
# chemin = "H:/Mon Drive/Recherche/Article 2 - Base de Donnees/7 - Data analysis/data/"

dcct = readxl::read_excel(paste0(chemin,
                                "cct_level_homogeneized_scored.xlsx")) 
dma = readxl::read_excel(paste0(chemin,
                                "ma_level_homogeneized_scored.xlsx"))
dcct$measure[dcct$measure == "SMC"] = "SMD"
dcct$measure[dcct$measure == "MC"] <- "MD"

dcct = dcct %>% filter(intervention_type == "Complementary")
dcct$factor = with(dcct, paste0(paper, "_", intervention_general, "_", outcome_general))
dma$Factor = with(dma, paste0(paper, "_", intervention_general, "_", outcome_general))

dat_error = view.errors.umbrella(dcct)
dcct[dat_error[dat_error$column_type_errors == "ERROR", "row_index"], "multiple_es"] <- "outcomes"
dat_error = view.errors.umbrella(dcct)
# rio::export(dcct, paste0("C:/Users/Corentin Gosling/drive_gmail/Recherche/Article 2 - Base de Donnees/7 - Data analysis/data/", "dat-CAM-analysis.txt"))

S1. PRIOR checklist

dat_checklist = readxl::read_excel()

S2. Deviations

S3. Search strategies

S4. Classification of age groups

  1. Homogeneous Preschool (<6yo): The 85th percentile of the distribution of the age is <6yo.

  2. Homogeneous School-age (6-12yo): The 25th percentile of the distribution of the age is >=6yo and the 75th percentile of the distribution of the age is < 13yo.

  3. Homogeneous Adolescents (13-19yo): The 25th percentile of the distribution of the age is >=13yo and the 75th percentile of the distribution of the age is < 20yo.

  4. Homogeneous Adults (>=20yo): The 15th percentile of the distribution of the age is >=18yo (because 18 is considered as the adult age in many countries).

S5. Included studies

S6. Excluded studies

SX. Primary analysis

SX-1. Data analysis

dcct_rr = dcct %>%
  group_by(factor) %>%
  filter(all(measure %in% c("OR", "RR")))

dcct_smd = dcct %>%
  group_by(factor) %>%
  filter(!all(measure %in% c("OR", "RR")))

if (
  !(all(dcct$factor %in% append(dcct_rr$factor, dcct_smd$factor)) &
  !(any(dcct_rr$factor %in% dcct_smd$factor) | any(dcct_smd$factor %in% dcct_rr$factor)))
) { stop("problem merging SMD and RR files")}


res_list_rr = umbrella(dcct_rr,
   verbose = FALSE,
   mult.level = TRUE, method.var = "PM", 
   r = 0.8, pre_post_cor = 0.5)

res_list_smd = umbrella(dcct_smd,
               verbose = FALSE,
               mult.level = TRUE, method.var = "REML", 
               r = 0.8, pre_post_cor = 0.5)
res_smd = summary(res_list_rr)
res_rr = summary(res_list_smd)
res_m_cred = bind_rows(res_smd, res_rr)
res_m_cred$Outcome = t(do.call(cbind, stringr::str_split(res_m_cred$Factor, "_")))[, 3]
res_m_cred$Intervention = t(do.call(cbind, stringr::str_split(res_m_cred$Factor, "_")))[, 2]
res_m_cred$"Meta-review" = t(do.call(cbind, stringr::str_split(res_m_cred$Factor, "_")))[, 1]
res_m = left_join(res_m_cred, dma, by="Factor") 
res_m$rob_num = paste0(round(res_m$rob), "%")
res_m$PICO_hom = grepl("Homo", res_m$PICO_amstar_precise, fixed=TRUE)

# View(res_m %>%
#        select(Factor, starts_with("available"), IN_meta, amstar.x, amstar.y, active_controls, mean_age))

SX-1. Certainty of evidence

Scoring components

RoB

res_m$down_rob = with(res_m,
    ifelse(rob > 75, 0, 
           ifelse(rob > 50, 1, 2)))
# View(res_m[, c("down_rob", "rob", "Factor")])

Inconsistency

# res_m$min_es = res_m$max_es = res_m$n_agreement = NA
# i = 0
# res_m$value = as.numeric(as.character(res_m$value))
# for (fact in unique(res_m$Factor)) {
#   i = i+1
#   # print(fact)
#   # print(i)
#   if (fact != res_m$Factor[i]) {
#     stop(paste0(fact, " != ", res_m$Factor[i], "(", i, ")"))
#   }
#   dat_sub = NA
#   dat_sub = res[[fact]]$x
#   values=pooled_es=NA
#   if (res_m$measure[i] == "G") {
#     values = dat_sub$value
#     pooled_es = res_m[res_m$Factor == fact, ]$value
#   } else {
#     values = log(dat_sub$value)
#     pooled_es = log(res_m[res_m$Factor == fact, ]$value)
#   }
#   res_m$min_es[i] = min(values)
#   res_m$max_es[i] = max(values)
#   res_m$n_agreement[i] = sum(sign(values) == sign(pooled_es)) / length(dat_sub$value)
# }
# res_m$down_het = with(res_m,
#     ifelse((I2 >= 50 & n_agreement < 0.90), 2,
#            ifelse(I2 >= 25 & n_agreement < 0.90, 1, 0)))

# View(res_m[, c("n_studies", "min_es", "max_es", "n_agreement", "I2", "down_het", "Factor")])

eq_range_or = c(0.8, 1.25)
eq_range_g = c(-0.1, 0.1)
meas_G = res_m$measure == "G"
meas_OR = res_m$measure %in% c("OR", "RR")

low_ci_pos = (meas_G & CI_lo_value > 0) | (meas_OR & CI_lo_value > 1)
low_ci_neg_range = (meas_G & CI_lo_value < 0) & (meas_G & CI_lo_value > eq_range_g[1]) |
                   (meas_OR & CI_lo_value < 1) & (meas_OR & CI_lo_value > eq_range_or[1])
low_ci_out_range = (meas_G & CI_lo_value < eq_range_g[1]) |
                   (meas_OR & CI_lo_value < eq_range_or[1])

low_pi_pos = (meas_G & PI_lo_value > 0) | (meas_OR & PI_lo_value > 1)
low_pi_neg_range = (meas_G & PI_lo_value < 0 & PI_lo_value > eq_range_g[1]) |
                   (meas_OR & PI_lo_value < 1) & (meas_OR & PI_lo_value > eq_range_or[1])
low_pi_out_range = (meas_G & PI_lo_value < eq_range_g[1]) |
                   (meas_OR & PI_lo_value < eq_range_or[1]) 

up_ci_neg = (meas_G & CI_up_value < 0) | (meas_OR & CI_up_value < 1)
up_ci_pos_range = (meas_G & CI_up_value > 0 & CI_up_value < eq_range_g[2]) |
                  (meas_OR & CI_up_value > 1) & (meas_OR & CI_up_value < eq_range_or[2])
up_ci_out_range = (meas_G & CI_up_value > eq_range_g[2]) |
                  (meas_OR & CI_up_value > eq_range_or[2])  

up_pi_neg = (meas_G & PI_up_value < 0) | (meas_OR & PI_up_value < 1)
up_pi_pos_range = (meas_G & PI_up_value > 0) & (meas_G & PI_up_value < eq_range_g[2]) |
                  (meas_OR & PI_up_value > 1) & (meas_OR & PI_up_value < eq_range_or[2])
up_pi_out_range = (meas_G & PI_up_value > eq_range_g[2]) |
                  (meas_OR & PI_up_value > eq_range_or[2]) 

# res_m$DOWN_HET = ifelse((low_ci_pos & low_pi_pos) | (up_ci_neg & up_pi_neg), 0,
#   ifelse((low_ci_neg_range & low_pi_neg_range) | (up_ci_pos_range & up_pi_pos_range), 0,
#   ifelse((low_ci_out_range & low_pi_out_range) | (up_ci_out_range & up_ci_out_range), 0,
#   ifelse((low_ci_pos & low_pi_neg_range) | (up_ci_neg & up_pi_pos_range), 1,
#   ifelse((low_ci_neg_range & low_pi_out_range) | (up_ci_pos_range & up_pi_out_range), 1,
#   ifelse((low_ci_pos & low_pi_out_range) | (up_ci_neg & up_ci_out_range), 2,
#   ifelse((low_ci_neg_range & low_pi_out_range) & (up_ci_pos_range & up_pi_out_range), 2, 999)))))))
res_m$down_het = 
  ifelse((low_ci_neg_range & low_pi_out_range) & (up_ci_pos_range & up_pi_out_range), 2, 
  ifelse((low_ci_pos & low_pi_out_range) | (up_ci_neg & up_ci_out_range), 2,
  ifelse((low_ci_neg_range & low_pi_out_range) | (up_ci_pos_range & up_pi_out_range), 1, 
  ifelse((low_ci_pos & low_pi_neg_range) | (up_ci_neg & up_pi_pos_range), 1,
  ifelse((low_ci_pos & low_pi_pos) | (up_ci_neg & up_pi_neg), 0,
  ifelse((low_ci_neg_range & low_pi_neg_range) | (up_ci_pos_range & up_pi_pos_range), 0,
  ifelse((low_ci_out_range & low_pi_out_range) | (up_ci_out_range & up_ci_out_range), 0,
  999)))))))

rows_smd = which(res_m$measure == "G")
rows_rr = which(res_m$measure != "G")

row_imp_pi = which(is.na(res_m$PI_lo_value) | is.na(res_m$PI_up_value))
res_m$down_het[row_imp_pi] <- 2


# View(cbind(res_m[, c("PI_lo_value", "CI_lo_value", "CI_up_value", "PI_up_value", 
#                "measure", "DOWN_HET", "n_studies")],
#                 low_ci_pos, low_ci_neg_range, low_ci_out_range,
#                low_pi_pos, low_pi_neg_range, low_pi_out_range,
#                up_ci_neg, up_ci_pos_range, up_ci_out_range,
#                up_pi_neg, up_pi_pos_range, up_pi_out_range)[rows_rr,])

Indirectness

res_m$down_ind =  with(res_m,
    ifelse(grepl("Mixed", res_m$PICO_amstar_ID, fixed = TRUE), 1,
    ifelse(res_m$available_control < 0.75 , 1, 0)
      # (active_controls >= 80 & active_controls <= 20) |
    ))

Imprecision

res_m$cross_low_high = (meas_G & (res_m$CI_lo_value < 0 & res_m$CI_up_value > 0.8) |
                          (res_m$CI_lo_value < -0.8 & res_m$CI_up_value > 0)) | 
  (meas_OR & (res_m$CI_lo_value < 1 & res_m$CI_up_value > 5) |
                          (res_m$CI_lo_value < 0.2 & res_m$CI_up_value > 1))

res_m$n_fail_detect_small_effects = 
  res_m$n_cases < 394 | res_m$n_controls < 394
res_m$n_fail_detect_large_effects = 
  res_m$n_cases < 64 | res_m$n_controls < 64

res_m$down_imp = ifelse(
  (res_m$cross_low_high & res_m$n_fail_detect_small_effects) |
    res_m$n_fail_detect_large_effects, 2,
  ifelse(res_m$cross_low_high | res_m$n_fail_detect_small_effects, 1, 0))

# View(res_m[, c( "down_imp", "ci_lo", "ci_up", "cross_low_high", "n_cases", "n_controls", "Factor")])

Publication bias

res_m$down_pubbias = ifelse(
  as.numeric(res_m$egger_p) < 0.10 |
  res_m$n_studies < 5 | 
  is.na(as.numeric(res_m$egger_p)), 1, 0)
## Warning in ifelse(as.numeric(res_m$egger_p) < 0.1 | res_m$n_studies < 5 | : NAs
## introduced by coercion

## Warning in ifelse(as.numeric(res_m$egger_p) < 0.1 | res_m$n_studies < 5 | : NAs
## introduced by coercion

Overal scoring

res_m$down_rob[is.na(res_m$down_rob)] <- 2
res_m$down_het[is.na(res_m$down_het)] <- 2
res_m$down_ind[is.na(res_m$down_ind)] <- 2
res_m$down_imp[is.na(res_m$down_imp)] <- 2
res_m$down_pubbias[is.na(res_m$down_pubbias)] <- 2

res_m$GRADE_num = rowSums(res_m[, 
        c("down_rob", "down_het", "down_ind", "down_imp", "down_pubbias")]) 
  
res_m$GRADE = with(res_m, ifelse(
  GRADE_num == 0, "High", 
  ifelse(GRADE_num %in% 1, "Moderate", 
         ifelse(GRADE_num %in% 2, "Low", 
                ifelse(GRADE_num >= 3, "Very low", NA)))))

res_m$Age = factor(res_m$Age, 
                   levels=c('< 6 yo', '6-12 yo', '13-19 yo', '>= 20 yo'))
res_m$GRADE = factor(res_m$GRADE, 
                   levels=c('High', 'Moderate', 'Low', 'Very low'))

res_m = res_m %>%
  arrange(Age, Intervention, GRADE, eG) 

# View(res_m %>%
#          filter(GRADE == "Moderate") %>% 
#        select(PICO_amstar_ID, IN_meta, n_studies, eG, eG_CI, 
#               down_rob, down_het, down_ind, down_imp, down_pubbias, 
#               I2, n_cases, n_controls, PI_eG))

colSums(res_m[, 
        c("down_rob", "down_het", "down_ind", 
          "down_imp", "down_pubbias")] == 2)
##     down_rob     down_het     down_ind     down_imp down_pubbias 
##           91           90            0          114            0
colSums(res_m[, 
        c("down_rob", "down_het", "down_ind", 
          "down_imp", "down_pubbias")] == 1)
##     down_rob     down_het     down_ind     down_imp down_pubbias 
##           58           37           80          127          166
res_m$paper = gsub("\\(child\\) ", "", res_m$paper)
res_m$paper = gsub("\\(adult\\) ", "", res_m$paper)
res_m$n_paper = length(unique(paste0(res_m$paper)))
interventions = sort(unique(res_m$intervention_general))
  
Mind_body_medicine = c("MUSIC", "SENS", "PHYS", "rTMS", "tCDS")

therapies_energetiques = c("ACUP");

Natural_Product_Based_Therapies  = c(
  "DIET", "FOLI", "HERB", "L-CARNIT", "L-CARNO", 
  "MELAT", "NAC", "OXYT", "PROB", "PUFA", 
  "SECRET", "SULFO", "VIT-D")
systemes_medicaux_alternatifs = c("AAI")

outcomes = sort(unique(res_m$outcome_general)); interventions; outcomes
##  [1] "AAI"      "ACUP"     "DIET"     "FOLI"     "HERB"     "L-CARNIT"
##  [7] "L-CARNO"  "MELAT"    "MUSIC"    "NAC"      "OXYT"     "PHYS"    
## [13] "PROB"     "PUFA"     "rTMS"     "SECRET"   "SENS"     "SULFO"   
## [19] "tDCS"     "VIT-D"
##  [1] "Acceptability"                   "Adaptative Functioning"         
##  [3] "ADHD symptoms"                   "Adverse events"                 
##  [5] "Anxiety"                         "Disruptive behaviors"           
##  [7] "Global cognition (IQ)"           "Language (Expressive skills)"   
##  [9] "Language (Overall skills)"       "Language (Receptive skills)"    
## [11] "Mood related symptoms"           "Overall ASD symptoms"           
## [13] "Quality of life"                 "Restricted/repetitive behaviors"
## [15] "Sensory Profile"                 "Sleep quality"                  
## [17] "Sleep quantity"                  "Social-communication"           
## [19] "Tolerability"

Primary analysis

SX-1. Synthesis - table

res_p = res_m %>% filter(IN_meta == 1)

# View(res_p %>%
#        group_by(PICO_amstar) %>%
#        summarise(n = n()))
DT::datatable(res_p,
      rownames = FALSE,
      extensions = 'Buttons',
      class = "display",
      options = list(
        # dom = c('t'),
        scrollX = TRUE,
        scrollCollapse = TRUE,

        scrollY = "500px",
        pageLength = nrow(res_p),
        columnDefs = list(
          list(width = '100px',
               targets = "_all"),
          list(className = 'dt-center',
               targets = "_all")),
        dom = c('tB'),
        buttons = c('copy', 'csv', 'excel','pdf')
        ))
# res_p %>%
#         filter(GRADE %in% c("Low", "Moderate") ) %>%
#         select(PICO_amstar, GRADE, Age, eG, p_value) %>%
#         arrange(Age, GRADE)
# res_m %>%
#         filter(Intervention == "OXYT" & Outcome == "Restricted/repetitive behaviors") %>%
#         select(PICO_amstar, GRADE, Age, eG, p_value, paper) %>%
#         arrange(Age, GRADE)

SX-1. Synthesis - Plots

Fig 2. — Low or Moderate GRADE rating

dat_low_or_higher = res_p %>% filter(GRADE %in% c("Low", "Moderate")) %>%
  mutate(age_pre = case_when(Age == "< 6 yo" ~ "Pre-school (<6 years old)",
                             Age == "6-12 yo" ~ "School-age (6-12 years old)",
                             Age == "13-19 yo" ~ "Adolescents (13-19 years old)",
                             Age == ">= 20 yo" ~ "Adults (>19 years old)")) %>%
  arrange(Age, GRADE, eG)
metaumbrella::forest(dat_low_or_higher,
       layout = "RevMan5",
       subgroup = "age_pre",
       subgroup.name = "",
       leftcols = c("Intervention", "Outcome", "GRADE", 
                    "n_studies",  "rob_num",  "I2", "effect.ci"),
       leftlabs = c("Intervention", "Outcome", "GRADE",
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Low or Moderate GRADE"
)

Core symptoms

res_core = res_p %>% filter(Outcome %in% ASD_symptoms)
metaumbrella::forest(res_core,
       layout = "RevMan5",
       subgroup = "Outcome",
       subgroup.name = "",
       leftcols = c("Intervention", "Age", "GRADE", 
                    "n_studies",  "rob_num",  "I2", "effect.ci"),
       leftlabs = c("Intervention", "Age", "GRADE",
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       # sortvar = res_p$Age,
       smlab = "Core ASD symptoms"
)

Safety

Fig 3. — Safety outcomes

forest(res_p %>% filter(Outcome %in% safety) %>%
  arrange(Age, GRADE, eG),
       title = "equivalent Risk Ratio (eRR)",
       xlim = c(0.1, 10),
              squaresize = 0.7,
       weight.study = "same",

       layout = "RevMan5",
       subgroup = "Outcome",
       xlab = "Equivalent Risk Ratio (eRR)",
       subgroup.name = "",
       leftcols = c("Intervention", "Age", "GRADE", 
                    "n_studies",  "rob_num",  "I2", "effect.ci"),
       leftlabs = c("Intervention", "Age", "GRADE", 
                    "n-studies", "Low\nRoB",   "I²", "eRR + 95% CI"),
       # sortvar = res_p$Age,
       smlab = "Safety"
)

Functioning

forest(res_p %>% filter(Outcome %in% functioning),
       layout = "RevMan5",
       subgroup = "Outcome",
       subgroup.name = "",
       leftcols = c("Intervention", "Age", "GRADE", 
                    "n_studies",  "rob_num",  "I2", "effect.ci"),
       leftlabs = c("Intervention", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       # sortvar = res_p$Age,
       smlab = "Daily functioning"
)

Comorbidities

forest(res_p %>% filter(Outcome %in% comorbidities),
       layout = "RevMan5",
       subgroup = "Outcome",
       subgroup.name = "",
       leftcols = c("Intervention", "Age", "GRADE", 
                    "n_studies",  "rob_num",  "I2", "effect.ci"),
       leftlabs = c("Intervention", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       # sortvar = res_p$Age,
       smlab = "Psychiatric comorbidities"
)

SX2. Overlapping

Overall picture

# Function to calculate interval overlap
interval_overlap_percentage <- function(lower1, upper1, lower2, upper2) {
  overlap = max(0, min(upper1, upper2) - max(lower1, lower2))
  total_length = (upper1 - lower1) + (upper2 - lower2) - overlap
  percentage_overlap = (overlap / total_length) * 100
  return(percentage_overlap)
}


res_over <- res_m %>%
  group_by(PICO_amstar) %>%
  filter(n() > 1) %>%
  mutate(n_SRMA = length(unique(paper))) %>%
  group_by(PICO_amstar) %>%
  mutate(
   eG_meta_IN = eG[IN_meta == 1],
   GRADE_meta_IN = GRADE[IN_meta == 1],
   p_value_meta_IN = p_value[IN_meta == 1],
    n = n(),
    n_SRMA = unique(n_SRMA),
    paper = collapsunique(paper),
    min_es = min(eG),
    max_es = max(eG),
    dif_es = max(eG) - min(eG),
    max_grade = head(sort(GRADE), 1),
    min_grade = tail(sort(GRADE), 1),
    prop_sig = sum(as.numeric(p_value) < 0.05) / n(),
    avg_percentage_overlap = mean(map_dbl(combn(n(), 2, simplify = FALSE), function(idx) {
    interval_overlap_percentage(
      ci_lo_g[idx[1]], ci_up_g[idx[1]], 
      ci_lo_g[idx[2]], ci_up_g[idx[2]])
  })),
    min_percentage_overlap = min(map_dbl(combn(n(), 2, simplify = FALSE), function(idx) {
    interval_overlap_percentage(
      ci_lo_g[idx[1]], ci_up_g[idx[1]], 
      ci_lo_g[idx[2]], ci_up_g[idx[2]])
  }))
  ) %>%
  arrange(PICO_amstar)
overlap_factors_disc = res_over %>% 
  # filter((PICO_hom & Age_precise == "Homogeneous") | !PICO_hom | IN_meta == 1) %>%
  filter(#abs(as.numeric(GRADE_meta_IN) - as.numeric(GRADE)) >= 2 | 
        (abs(eG_meta_IN - eG) >= 0.30 & 
           sum(as.numeric(p_value_meta_IN) < 0.05,
               as.numeric(p_value) < 0.05) == 1) |
          as.numeric(min_percentage_overlap) < 10 |
          as.numeric(avg_percentage_overlap) < 20) %>%
  select(PICO_amstar, Factor)

res_over_disc = res_over %>%
  filter(PICO_amstar %in% overlap_factors_disc$PICO_amstar)


res_test = res_m %>%
 group_by(PICO_amstar) %>%
  filter(n() > 1) %>%
  mutate(PICO_hom = grepl("Homo", PICO_amstar_precise, fixed=TRUE)) %>%
  filter((PICO_hom & Age_precise == "Homogeneous") | !PICO_hom | IN_meta == 1) %>%
  group_by(PICO_amstar) %>%
  filter(n() > 1) %>%
  mutate(eG_meta_IN = eG[IN_meta == 1],
         GRADE_meta_IN = GRADE[IN_meta == 1],
         p_value_meta_IN = p_value[IN_meta == 1]) %>%
  ungroup() %>% rowwise() %>%
  filter(
    (abs(eG - eG_meta_IN) >= 0.30 & as.numeric(p_value) != as.numeric(p_value_meta_IN)) | 
    (GRADE == GRADE_meta_IN & as.numeric(p_value) != as.numeric(p_value_meta_IN)) |
    abs(as.numeric(GRADE) - as.numeric(GRADE_meta_IN)) >= 2 |
           IN_meta == 1) %>%
  # select(c(1:15, PICO_amstar, PICO_amstar_precise, IN_meta, amstar.x, Age, Age_precise, eG_meta_IN, GRADE, GRADE_meta_IN)) %>%
  arrange(PICO_amstar)
         
         



propC = function(x, R=0) {
  x_non_na = x[!is.na(x)]
  round(sum(x_non_na)/length(x) * 100, R)
}

DT::datatable(res_over)
# identify, for each 

Fig 4. Overlap

res_m$inter_out = paste0(res_m$Intervention, " - ", res_m$Outcome)
metaumbrella::forest(res_m %>% filter(Factor %in% res_over_disc$Factor),
       layout = "RevMan5",
       squaresize = 0.7,
       weight.study = "same",
       subgroup = "inter_out",
       subgroup.name = "",
       leftcols = c("paper", "Age", "GRADE", 
                    "n_studies",  "rob_num", "I2", "effect.ci"),
       leftlabs = c("Paper", "Age", "GRADE",
                    "n-studies", "Low RoB", "I²", "eSMD + 95% CI"),
       # sortvar = res_p$Age,
       smlab = "Key overlapping situations"
)

Mind-Body Medicine

MUSIC

forest(res_m %>% filter(Intervention %in% "MUSIC"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR (key items)", "AMSTAR (Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (MUSIC)")

SENS

forest(res_m %>% filter(Intervention %in% "SENS"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR (key items)", "AMSTAR (Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (SENS)")

PHYS

forest(res_m %>% filter(Intervention %in% "PHYS"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR (key items)", "AMSTAR (Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (PHYS)")

rTMS

forest(res_m %>% filter(Intervention %in% "rTMS"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR (key items)", "AMSTAR (Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (rTMS)")

tCDS

forest(res_m %>% filter(Intervention %in% "tDCS"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (tDCS)")

Energy Medicine

ACUP

forest(res_m %>% filter(Intervention %in% "ACUP"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (ACUP)")

Natural Product

Hormones

MELAT
forest(res_m %>% filter(Intervention %in% "MELAT"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (MELAT)")

SECRET
forest(res_m %>% filter(Intervention %in% "SECRET"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (SECRET)")

OXYT
forest(res_m %>% filter(Intervention %in% "SECRET"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (SECRET)")

Supplementation

Natural_Product_Based_Therapies = c( ““,”“,”“,”“,”“,”“,”“,”“,”“,”“,”“,”“,”“,”“)

FOLI
forest(res_m %>% filter(Intervention %in% "FOLI"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (FOLI)")

HERB
forest(res_m %>% filter(Intervention %in% "HERB"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (HERB)")

L-CARNIT
forest(res_m %>% filter(Intervention %in% "L-CARNIT"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (L-CARNIT)")

L-CARNO
forest(res_m %>% filter(Intervention %in% "L-CARNO"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (L-CARNO)")

NAC
forest(res_m %>% filter(Intervention %in% "NAC"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (NAC)")

PROB
forest(res_m %>% filter(Intervention %in% "PROB"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (L-CARNIT)")

SULFO
forest(res_m %>% filter(Intervention %in% "SULFO"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (SULFO)")

PUFA
res_pufa = res_m %>% filter(Intervention %in% "PUFA")

forest(res_pufa %>% filter(!outcome_general %in% safety),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (PUFA)")

forest(res_pufa %>% filter(outcome_general %in% safety),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (PUFA)")

VIT-D
forest(res_m %>% filter(Intervention %in% "VIT-D"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (VIT-D)")

Diets

forest(res_m %>% filter(Intervention %in% "DIET"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (DIET)")

Whole Medical Systems

AAI
forest(res_m %>% filter(Intervention %in% "AAI"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (AAI)")

SX. Writing

Characteristics

paste0("The ", unique(res_m$n_paper), " retained reports explored the effects of ", 
       length(unique(res_m$intervention_general)), " intervention types, on at least one of the ", 
       length(unique(res_m$outcome_general)), " possible outcomes.",
       " These reports generated a total of ", nrow(res_m), " meta-analyses",
       " that synthesized the evidence from more than ",
       length(unique(paste0(dcct$author, dcct$year))), " CCTs, and that explored ", length(unique(res_m$PICO_amstar)), " unique PICOs.")
## [1] "The 52 retained reports explored the effects of 20 intervention types, on at least one of the 19 possible outcomes. These reports generated a total of 242 meta-analyses that synthesized the evidence from more than 202 CCTs, and that explored 142 unique PICOs."
res_year = res_m %>%
  group_by(paper) %>%
  slice(1)
paste0("These reports were, in their vast majority (", 
       propC(res_year$year_meta_num >= 2018), "%) published after 2018.")
## [1] "These reports were, in their vast majority (87%) published after 2018."
dat_pico_overlap = res_m %>%
  group_by(PICO_amstar, Intervention, Outcome, Age) %>%
  summarise(n=n()) %>%
  filter(n > 1)
## `summarise()` has grouped output by 'PICO_amstar', 'Intervention', 'Outcome'.
## You can override using the `.groups` argument.
paste0("A total of ", nrow(dat_pico_overlap), " PICOs (", round(nrow(dat_pico_overlap)/length(unique(res_m$PICO_amstar)) * 100), "%) was explored in several overlapping meta-analyses. The PICOs with the highest number of overlapping meta-analyses were investigating the effects of ", 
        unique(dat_pico_overlap[dat_pico_overlap$n >= 6, ]$Intervention), " on ", tolower(paste(paste0(paste0(unique(dat_pico_overlap[dat_pico_overlap$n >= 6, ]$Outcome)), " (n=",
paste0(unique(dat_pico_overlap[dat_pico_overlap$n >= 6, ]$n)), ")")
, collapse=", ")))
## [1] "A total of 52 PICOs (37%) was explored in several overlapping meta-analyses. The PICOs with the highest number of overlapping meta-analyses were investigating the effects of PUFA on adhd symptoms (n=7), disruptive behaviors (n=8), restricted/repetitive behaviors (n=7), social-communication (n=8)"
# paste0("The mean number of overlapping meta-analyses on the same PICO was ",
#       round(mean(dat_pico$n), 2), " (minimum = ", 
#       min(dat_pico$n), ", maximum = ",
#       max(dat_pico$n), ").")

# paste0(length(unique(res_m$intervention_general)), " intervention types (", sort(collapsunique(res_m$intervention_general)), ")")
# paste0(length(unique(res_m$outcome_general)), " outcomes (", sort(collapsunique(res_m$outcome_general)), ")")

Description of the meta-analyses

res_interventions = res_p %>%
  group_by(Intervention) %>%
  select(Outcome) %>%
  distinct() %>%
  summarise(n = n(),
            outcome = collapsunique(Outcome)) %>%
  arrange(desc(n))
## Adding missing grouping variables: `Intervention`
nrowInter = nrow(res_interventions)

table(res_p$amstar_rank)
## 
## Critically low           High            Low 
##             32             23             87
dat_amstar = as.data.frame(table(res_p %>% group_by(paper) %>% slice(1) %>% ungroup%>% select(amstar_rank)))
dat_amstar_full = as.data.frame(table(res_p %>%  select(amstar_rank)))

paste0("Thus, ", length(unique(res_m$PICO_amstar)), " meta-analyses, derived from ", length(unique(res_p$paper)), " reports, were finally included in the primary analysis after discarding overlapping meta-analyses.")
## [1] "Thus, 142 meta-analyses, derived from 31 reports, were finally included in the primary analysis after discarding overlapping meta-analyses."
paste0("According to the AMSTAR-2 scoring, ",
        dat_amstar_full[dat_amstar_full$Var1=="High", "Freq"],
        " PICOs were evaluated by meta-analyses of high quality (",
       round(dat_amstar_full[dat_amstar_full$Var1=="High", "Freq"] / sum(dat_amstar_full$Freq) * 100), "%), ",
        dat_amstar_full[dat_amstar_full$Var1=="Low", "Freq"], 
        " of low quality (",
        round(dat_amstar_full[dat_amstar_full$Var1=="Low", "Freq"] / sum(dat_amstar_full$Freq) * 100), "%), and ", dat_amstar_full[dat_amstar_full$Var1=="Critically low", "Freq"],
        " of critically low quality (", 
        round(dat_amstar_full[dat_amstar_full$Var1=="Critically low", "Freq"] / sum(dat_amstar_full$Freq) * 100), "%), ",
       
        dat_amstar_full[dat_amstar_full$Var1=="Moderate", "Freq"], " moderate")
## [1] "According to the AMSTAR-2 scoring, 23 PICOs were evaluated by meta-analyses of high quality (16%), 87 of low quality (61%), and 32 of critically low quality (23%),  moderate"
paste0("According to the AMSTAR-2 scoring, ",
        dat_amstar[dat_amstar$Var1=="High", "Freq"],
        " of these ", length(unique(res_p$paper)), 
        " reports were of high quality, ",
        dat_amstar[dat_amstar$Var1=="Low", "Freq"], 
        " of low quality, and ",
        dat_amstar[dat_amstar$Var1=="Critically low", "Freq"],
        " of critically low quality. ",
        dat_amstar[dat_amstar$Var1=="Moderate", "Freq"], " moderate")
## [1] "According to the AMSTAR-2 scoring, 4 of these 31 reports were of high quality, 9 of low quality, and 18 of critically low quality.  moderate"
paste0("The interventions that covered the largest number of outcomes were ", paste(paste0(res_interventions$Intervention[1:5], " (n=", res_interventions$n[1:5], "/",length(unique(res_p$Outcome)), ")"), collapse=", "))
## [1] "The interventions that covered the largest number of outcomes were OXYT (n=11/19), PUFA (n=11/19), AAI (n=10/19), ACUP (n=9/19), NAC (n=9/19)"
# paste0("The interventions that covered the smallest number of outcomes were ", paste(paste0(res_interventions$Intervention[(nrowInter-8):nrowInter], " (n=", res_interventions$n[(nrowInter-8):nrowInter], "/",length(unique(res_p$Outcome)), ")"), collapse=", "))

res_outcomes = res_p %>%
  group_by(Outcome) %>%
  select(Intervention) %>%
  distinct() %>%
  summarise(n = n(),
            outcome = collapsunique(Intervention)) %>%
  arrange(desc(n))
## Adding missing grouping variables: `Outcome`
nrowOutcome = nrow(res_outcomes)
paste0("and the outcomes that were covered by the largest number of interventions all regarded the efficacy of the interventions on core ASD symptoms: ", paste(paste0(tolower(res_outcomes$Outcome[1:5]), " (n=", res_outcomes$n[1:5], "/",length(unique(res_m$Intervention)), ")"), collapse=", "))
## [1] "and the outcomes that were covered by the largest number of interventions all regarded the efficacy of the interventions on core ASD symptoms: social-communication (n=17/20), overall asd symptoms (n=16/20), restricted/repetitive behaviors (n=13/20), disruptive behaviors (n=11/20), adhd symptoms (n=9/20)"
res_adverse = res_p %>%
  filter(Outcome %in% safety) %>%
  group_by(Outcome, Intervention) %>%
  distinct() %>%
  summarise(n = n()) %>%
  group_by(Intervention) %>%
  summarise(n = n())
## `summarise()` has grouped output by 'Outcome'. You can override using the
## `.groups` argument.
paste0("Only ", nrow(res_adverse[res_adverse$n==3, ]),
       " interventions were tested for all our safety outcomes, namely acceptability, tolerability, and the presence of adverse effects (", paste(sort(unlist(res_adverse[res_adverse$n==3, "Intervention"])), collapse=", "), "). Some others were tested for one or two of these outcomes (", paste(sort(unlist(res_adverse[res_adverse$n!=3, "Intervention"])), collapse=", "), ").")
## [1] "Only 4 interventions were tested for all our safety outcomes, namely acceptability, tolerability, and the presence of adverse effects (MELAT, NAC, OXYT, PUFA). Some others were tested for one or two of these outcomes (L-CARNO, PROB, SULFO, VIT-D)."
# sort(unique(lapply(res_outcomes[res_outcomes$Outcome %in% safety, "outcome"], function(x) unlist(strsplit(x, ", ")))[[1]]))
# res_outcomes[res_outcomes$Outcome %in% safety, ]
# sort(unique(lapply(res_outcomes[res_outcomes$Outcome %in% safety, "outcome"], function(x) unlist(strsplit(x, ", ")))[[1]]))



# paste0("The outcomes that were covered by the least number of interventions were ", paste(paste0(res_outcomes$Outcome[(nrowOutcome-5):nrowOutcome], " (n=", res_outcomes$n[(nrowOutcome-5):nrowOutcome], "/",length(unique(res_p$Intervention)), ")"), collapse=", "))

paste0("The median number of RCTs per meta-analysis was ",
      quantile(res_p$n_studies, 0.5), " (IQR = [", 
      quantile(res_p$n_studies, 0.25), ", ",
      quantile(res_p$n_studies, 0.75), "]), and the median number of participants per meta-analysis was ",
      quantile(res_p$total_n, 0.5), " (IQR = [", 
      quantile(res_p$total_n, 0.25), ", ",
      quantile(res_p$total_n, 0.75), "])")
## [1] "The median number of RCTs per meta-analysis was 4 (IQR = [2, 5]), and the median number of participants per meta-analysis was 154 (IQR = [92, 236])"
paste0("The median percentage of female per meta-analysis was ",
      round(quantile(res_p$perc_female, 0.5, na.rm=TRUE)), "% (IQR = [", 
      round(quantile(res_p$perc_female, 0.25, na.rm=TRUE)), ", ",
      round(quantile(res_p$perc_female, 0.75, na.rm=TRUE)), "])")
## [1] "The median percentage of female per meta-analysis was 14% (IQR = [10, 19])"
tab_Age = as.data.frame(table(res_p$Age))
paste0("A total of ",
       tab_Age[tab_Age$Var1=="< 6 yo", "Freq"], 
       " meta-analyses (",
       round(tab_Age[tab_Age$Var1=="< 6 yo", "Freq"] / nrow(res_p) * 100),
       "%) involved very young children (<6 years old), ", 
        tab_Age[tab_Age$Var1== "6-12 yo", "Freq"], 
       " meta-analyses (",
      round(tab_Age[tab_Age$Var1 == "6-12 yo", "Freq"] / nrow(res_p) * 100), "%) involved school-aged children (6-12 years old), and ",
        tab_Age[tab_Age$Var1== "13-19 yo", "Freq"], " meta-analyses (",
round(tab_Age[tab_Age$Var1 == "13-19 yo", "Freq"] / nrow(res_p) * 100), "%) involved adolescents (13-19 years old), and ",
        tab_Age[tab_Age$Var1== ">= 20 yo", "Freq"], " meta-analyses (",
round(tab_Age[tab_Age$Var1 == ">= 20 yo", "Freq"] / nrow(res_p) * 100), "%) involved adults (>= 20 years old)")
## [1] "A total of 26 meta-analyses (18%) involved very young children (<6 years old), 85 meta-analyses (60%) involved school-aged children (6-12 years old), and 15 meta-analyses (11%) involved adolescents (13-19 years old), and 16 meta-analyses (11%) involved adults (>= 20 years old)"
tab_IQ = as.data.frame(table(res_p$IQ))
paste0("A total of ",
       tab_IQ[tab_IQ$Var1=="Average (80-119)", "Freq"], " meta-analyses (",
       round(tab_IQ[tab_IQ$Var1=="Average (80-119)", "Freq"] / nrow(res_p) * 100), "%) focused on participants with overall cognitive functioning within the norm (IQ = [80-119]), ", 
        sum(tab_IQ[tab_IQ$Var1 %in% c("Limit (70-79)", "Low (< 70)"), "Freq"]), " meta-analyses (",
round(sum(tab_IQ[tab_IQ$Var1 %in% c("Limit (70-79)", "Low (< 70)"), "Freq"]) / nrow(res_p) * 100), "%) focused on participant with moderate or important cognitive difficulties (IQ < 80), and ",
        nrow(res_p) - sum(tab_IQ$Freq), " (",
       round(100-sum(tab_IQ$Freq) / nrow(res_p) * 100), "%) did not provide information on the cognitive functioning of the participants.")
## [1] "A total of 62 meta-analyses (44%) focused on participants with overall cognitive functioning within the norm (IQ = [80-119]), 22 meta-analyses (15%) focused on participant with moderate or important cognitive difficulties (IQ < 80), and 58 (41%) did not provide information on the cognitive functioning of the participants."
res_large = res_p %>% filter(eG > 0.5)
res_sig = res_p %>% filter(as.numeric(p_value) < 0.05)
res_effect = res_p %>% filter(eG > 0 & (eG > 0.5 | as.numeric(p_value) < 0.05))

Effect size magnitude

paste0("The overall results indicated that ", 
       sum(res_p$eG > 0), " meta-analyses (", propC(res_p$eG > 0), 
       "%) had a pooled effect size in favor of the experimental group (SMD > 0, OR/RR > 1). More precisely, ", 
       nrow(res_large), " (",
       round(nrow(res_large) / nrow(res_p) * 100) , "%) had a magnitude that can be considered as moderate or large (SMD > 0.50), and ",
       nrow(res_sig), " (",
       round(nrow(res_sig) / nrow(res_p) * 100) , "%) had a statistically significant p-value (p < 0.05); (",
       round(nrow(res_effect) / nrow(res_p) * 100), "%")
## [1] "The overall results indicated that 109 meta-analyses (77%) had a pooled effect size in favor of the experimental group (SMD > 0, OR/RR > 1). More precisely, 26 (18%) had a magnitude that can be considered as moderate or large (SMD > 0.50), and 35 (25%) had a statistically significant p-value (p < 0.05); (28%"
paste0("It should be acknowledged that the ",
       nrow(res_effect), " meta-analyses with moderate to large pooled effect sizes and/or statistically significant results are not without limitations, given that only ",
       sum(res_effect$rob > 75), " of them (", propC(res_effect$rob > 75), "%) had more than 75% of their participants coming from low-risk studies.",
       " Moreover, ",
       sum(res_effect$I2 >= 25), " (", propC(res_effect$I2 >= 25),
       "%) of these meta-analyses had a moderate or large inconsistency (I² >= 25%) and only ",
       sum(as.numeric(res_effect$PI_lo_g) > 0 | as.numeric(res_effect$PI_up_g) < 0, na.rm=TRUE), " (", 
       propC(as.numeric(res_effect$PI_lo_g) > 0 | as.numeric(res_effect$PI_up_g) < 0), "%) of the meta-analyses had a 95% prediction interval that excluded the null.")
## [1] "It should be acknowledged that the 40 meta-analyses with moderate to large pooled effect sizes and/or statistically significant results are not without limitations, given that only 7 of them (18%) had more than 75% of their participants coming from low-risk studies. Moreover, 18 (45%) of these meta-analyses had a moderate or large inconsistency (I² >= 25%) and only 4 (10%) of the meta-analyses had a 95% prediction interval that excluded the null."
# . Moreover, in ", propC(res_effect$n_agreement < 0.8), "% of the meta-analyses, more than 25% of the individual effect sizes were of the opposite sign to the pooled effect size

intervention_ns = unique(as.character(unlist(res_p %>%
  filter(!Outcome %in% safety) %>%
  group_by(Intervention) %>%
  filter(all(as.numeric(p_value) >= 0.05)) %>%
  select(Intervention))))

intervention_s = unique(as.character(unlist(res_p %>%
  filter(!Outcome %in% safety) %>%
  group_by(Intervention) %>%
  filter(any(as.numeric(p_value) < 0.05)) %>%
  select(Intervention))))

if (length(unique(append(intervention_ns, intervention_s))) != 20) {
   stop("problem identifying ns interventions")
}
paste0("A total of ", length(unique(intervention_ns)), " interventions presented no evidence of efficacy on any of the target outcomes (", paste(sort(intervention_ns), collapse=", "), ").")
## [1] "A total of 7 interventions presented no evidence of efficacy on any of the target outcomes (FOLI, L-CARNIT, L-CARNO, PROB, SECRET, SULFO, VIT-D)."
# View(res_p %>%
#        filter(Outcome == "Overall ASD symptoms"))
       
paste0(
  ifelse(nrow(res_p %>% filter(eG < 0  & as.numeric(p_value) < 0.05)), " there is a detrimental effect", "there is no detrimental effect" 
))
## [1] "there is no detrimental effect"

GRADING

nrow(res_p[res_p$GRADE == "High", ])
## [1] 0
dat_grade = as.data.frame(table(res_p$GRADE))
paste0("The assessment of the certainty of evidence allowed the emergence of a very clear pattern. Of the ", nrow(res_p), " unique PICOs retained in our main analysis, ", 
       dat_grade[dat_grade$Var1 == "High", "Freq"],
       " reached a 'High' GRADE rating, ", 
       dat_grade[dat_grade$Var1 == "Moderate", "Freq"],
       " reached a 'Moderate' rating, ",
       dat_grade[dat_grade$Var1 == "Low", "Freq"],
       " (", round(dat_grade[dat_grade$Var1 == "Low", "Freq"] / sum(dat_grade$Freq) * 100),"%) reached a 'Low' rating, and ",
       dat_grade[dat_grade$Var1 == "Very low", "Freq"],
       " (", round(dat_grade[dat_grade$Var1 == "Very low", "Freq"] / sum(dat_grade$Freq) * 100), "%) reached a 'Very low' rating. ")
## [1] "The assessment of the certainty of evidence allowed the emergence of a very clear pattern. Of the 142 unique PICOs retained in our main analysis, 0 reached a 'High' GRADE rating, 5 reached a 'Moderate' rating, 14 (10%) reached a 'Low' rating, and 123 (87%) reached a 'Very low' rating. "
res_mod = res_p %>% filter(GRADE == "Moderate")
paste0("The ", nrow(res_mod), " PICOs with a 'Moderate' GRADE rating explored the effects of ", paste(paste0(res_mod$intervention_general, " on ", res_mod$outcome_general, " in ", res_mod$Age, " (SMD=", res_mod$eG, ", 95% CI=",
                                                                                                             res_mod$eG_CI, ")"), collapse= ", "))
## [1] "The 5 PICOs with a 'Moderate' GRADE rating explored the effects of OXYT on Overall ASD symptoms in 6-12 yo (SMD=-0.016, 95% CI=[-0.209, 0.178]), OXYT on Acceptability in >= 20 yo (SMD=-0.025, 95% CI=[-0.561, 0.511]), OXYT on Overall ASD symptoms in >= 20 yo (SMD=0.008, 95% CI=[-0.269, 0.284]), OXYT on Tolerability in >= 20 yo (SMD=0.141, 95% CI=[-0.549, 0.831]), OXYT on Restricted/repetitive behaviors in >= 20 yo (SMD=0.395, 95% CI=[0.137, 0.653])"
low_mod_sig = dat_low_or_higher[dat_low_or_higher$eG > 0 & as.numeric(dat_low_or_higher$p_value) < 0.05, ]
paste0(" ", nrow(low_mod_sig))
## [1] " 2"
paste0("More precisely, when considering only PICOs with a ‘Low’ or ‘Moderate’ rating, only ", paste(paste0(low_mod_sig$intervention_general, " on ", low_mod_sig$outcome_general, " in ", low_mod_sig$Age, " (SMD=", low_mod_sig$eG, ", 95% CI=",
                                                                                                             low_mod_sig$eG_CI, ")"), collapse= ", "))
## [1] "More precisely, when considering only PICOs with a ‘Low’ or ‘Moderate’ rating, only MUSIC on Overall ASD symptoms in < 6 yo (SMD=0.211, 95% CI=[0.023, 0.399]), OXYT on Restricted/repetitive behaviors in >= 20 yo (SMD=0.395, 95% CI=[0.137, 0.653])"

Overlap

paste0("Among the ", length(unique(res_over$PICO_amstar)), " PICOs that were explored by at least 2 overlapping meta-analyses, the median overlap of the 95% CI of the pooled effect sizes was ",
       round(quantile(res_over$avg_percentage_overlap, 0.5)),
       "% (IQR=[", 
       round(quantile(res_over$avg_percentage_overlap, 0.25)),
       "%, ",
       round(quantile(res_over$avg_percentage_overlap, 0.75)),
       "%]). Only ", 
       sum(res_over$min_grade == res_over$max_grade), " (", 
       round(sum(res_over$min_grade == res_over$max_grade) / nrow(res_over) * 100), 
       "%) of PICOs consistently achieved a similar GRADE rating and only ", 
       sum(res_over$prop_sig %in% c(0,1)), " (", 
       round(sum(res_over$prop_sig %in% c(0,1)) / nrow(res_over) * 100), 
       "%) of PICOs consistently achieved a similar statistical significance status.") 
## [1] "Among the 52 PICOs that were explored by at least 2 overlapping meta-analyses, the median overlap of the 95% CI of the pooled effect sizes was 52% (IQR=[43%, 65%]). Only 82 (54%) of PICOs consistently achieved a similar GRADE rating and only 97 (64%) of PICOs consistently achieved a similar statistical significance status."
fact_sig = res_p %>% filter(as.numeric(p_value) < 0.05) %>% select(PICO_amstar)

res_sig_overlap = res_over %>% filter(PICO_amstar %in% as.character(unlist(fact_sig)))

paste0("Among the ", nrow(res_sig_overlap), " PICOs with statistically significant results in our primary analysis and at least one overlapping meta-analysis, we found that ",
       nrow(res_sig_overlap %>% filter(!prop_sig %in% c(0, 1))),
       " (",
       round(nrow(res_sig_overlap %>% filter(!prop_sig %in% c(0, 1))) / nrow(res_sig_overlap) * 100), "%) presented at least one overlapping meta-analysis with non-statistically significant results")
## [1] "Among the 49 PICOs with statistically significant results in our primary analysis and at least one overlapping meta-analysis, we found that 39 (80%) presented at least one overlapping meta-analysis with non-statistically significant results"
paste0("We found ", length(unique(res_over_disc$PICO_amstar)), " PICOs that presented a large discrepancy regarding either the magnitude of the estimated effect (i.e., at least two pooled effect sizes varying by more than SMD >= 0.30 with different statistical significance status), or regarding the confidence or GRADE status differing by two points (e.g., two meta-analyses with 'Very low' versus 'Moderate' GRADE rating; Fig 4.).")
## [1] "We found 6 PICOs that presented a large discrepancy regarding either the magnitude of the estimated effect (i.e., at least two pooled effect sizes varying by more than SMD >= 0.30 with different statistical significance status), or regarding the confidence or GRADE status differing by two points (e.g., two meta-analyses with 'Very low' versus 'Moderate' GRADE rating; Fig 4.)."

OXYT restrict

dat_oxyt_rest = dcct %>% 
   filter(intervention_general == "NAC" &
          first_author_meta %in% c("Iffland", "Siafis (child)") &
          outcome_general == "Disruptive behaviors") %>%
  arrange(author, year) %>%
  mutate(TE = ifelse(!is.na(reverse_es), -as.numeric(as.character(value)), value),
         TE_lo = ifelse(!is.na(reverse_es), -as.numeric(as.character(ci_up)), ci_lo),
         TE_up = ifelse(!is.na(reverse_es), -as.numeric(as.character(ci_lo)), ci_up), 
         ID = paste0(author, " (", year, ")"))

met_oxyt_res = meta::metagen(TE = as.numeric(TE),
   lower = as.numeric(TE_lo),
   upper = as.numeric(TE_up),
   subgroup = ID,
   subgroup.name = "Study",
   studlab = meta_review,
   data = dat_oxyt_rest)


m2 <- meta::rob(RoB_SeqGen, RoB_Conceal, RoB_BlindPers, RoB_BlindOut, RoB_Attrition, RoB_Reporting, overall = rob, data = met_oxyt_res, tool = "rob1")

meta::forest(m2,
       layout = "RevMan5", common = FALSE, random = FALSE, 
       leftcols = c("studlab", "effect.ci"),
       leftlabs = c("Study","SMD [95% CI]"),
       col.square = "gray", 
  col.study = "black", col.square.lines = "black",
  overall = FALSE, hetstat = FALSE)

#        subgroup = "author",
#        subgroup.name = ""
# metaumbrella::forest(focus_oxyt_restrict_adult,
#        layout = "RevMan5",
#        subgroup = "inter_out",
#        subgroup.name = "",
#        leftcols = c("paper", "Age", "Age_precise", "GRADE", 
#                     "n_studies",  "rob", "active_controls", "I2", "effect.ci"),
#        leftlabs = c("Paper", "Age", "Age_precise", "GRADE",
#                     "n-studies", "Low RoB", "Active controls",  "I²", "eSMD + 95% CI"),
#        # sortvar = res_p$Age,
#        smlab = "Key overlapping situations"
# )
res_over_sig <- res_m %>%
  group_by(PICO_amstar) %>%
  filter(n() > 1) %>%
  mutate(eG_meta_IN = eG[IN_meta == 1],
         GRADE_meta_IN = GRADE[IN_meta == 1],
         p_value_meta_IN = p_value[IN_meta == 1]) %>%
  filter(as.numeric(p_value_meta_IN) < 0.05) %>%
  summarise(
    n = n(),
    paper = collapsunique(paper),
    factor = collapsunique(Factor),
    min_es = min(eG),
    max_es = max(eG),
    dif_es = max(eG) - min(eG),
    max_grade = head(sort(GRADE), 1),
    min_grade = tail(sort(GRADE), 1),
    prop_sig = sum(as.numeric(p_value) < 0.05) / n(),
    avg_percentage_overlap = mean(map_dbl(combn(n(), 2, simplify = FALSE), function(idx) {
    interval_overlap_percentage(
      ci_lo_g[idx[1]], ci_up_g[idx[1]], 
      ci_lo_g[idx[2]], ci_up_g[idx[2]])
  })))


nrow(res_over_sig %>% filter(prop_sig < 1)) / nrow(res_over_sig)
## [1] 0.6875
res_discrepant = res_m %>% filter(PICO_amstar %in% res_over_disc$PICO_amstar)
metaumbrella::forest(res_discrepant,
       layout = "RevMan5",
       subgroup = "PICO_amstar",
       subgroup.name = "",
       leftcols = c("Intervention", "Outcome", "Age_precise", "paper", 
                    "GRADE", 
                    "n_studies",  "rob_num",  
                    "I2", "effect.ci"),
       leftlabs = c("Intervention", "Outcome", "Homogeneous Age", "Paper", 
                    "GRADE",
                    "n-studies", "Low RoB",   
                    "I²", "eSMD + 95% CI"),
       # sortvar = res_p$Age,
       smlab = "Discrepant PICOs"
)

propC((res_over$min_grade == "Very low" & res_over$max_grade == "Moderate") |
       res_over$dif_es>= 0.50 & !res_over$prop_sig %in% c(0, 1))
## [1] 29
paste0("The median overlap per meta-analysis was ",
      quantile(res_over$avg_percentage_overlap, 0.5), " (IQR = [", 
      quantile(res_over$avg_percentage_overlap, 0.25), ", ",
      quantile(res_over$avg_percentage_overlap, 0.75), "])")
## [1] "The median overlap per meta-analysis was 52.3123727984725 (IQR = [43.3679354094579, 65.2751482834739])"
res_m_cred = metaumbrella::add.evidence(res,
                           criteria = "Personalized",
                           class_I = c(total_n = 500, p_value = 1e-3, I2 = 50,
                                       egger_p = .05, pi = "notnull",
                                       rob = 75),
                           class_II = c(total_n = 350, p_value = 1e-3,
                                        largest_CI = "notnull",
                                        rob = 50),
                           class_III = c(total_n = 200, p_value = 1e-3),
                           class_IV = c(p_value = 5e-2), verbose = FALSE)